function pscb = meanpsctrans(b, mb, mt)
% meanpsctrans  - perform PSC transformation after regression
%
% FORMAT:       pscb = meanpsctrans(b, mb [, mt])
%
% Input fields:
%
%       b           beta images
%       mb          mean confound beta image
%       mt          mask threshold (defaut: 0.25)
%
% Output fields:
%
%       pscb        PSC transformed beta image (masked)

% Version:  v0.7g
% Build:    9040222
% Date:     Apr-02 2009, 10:18 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
   ~isnumeric(b) || ...
   ~isnumeric(mb)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument.' ...
    );
end
if nargin < 3 || ...
   ~isa(mt, 'double') || ...
    numel(mt) ~= 1 || ...
    isinf(mt) || ...
    isnan(mt) || ...
    mt < 0 || ...
    mt > 0.5
    mt = 0.25;
end
sb = size(b);
if prod(sb(1:end-1)) ~= numel(mb)
    error( ...
        'BVQXtools:BadArgument', ...
        'Last dimension of b must be image dimension.' ...
    );
end

% reshape mean beta to make sure and get repmat argument
if numel(sb) == 2
    mb = reshape(mb, sb(1), 1);
    rma = [1, sb(2)];
else
    mb = reshape(mb, sb(1:end-1));
    rma = [ones(1, numel(sb) - 1), sb(end)];
end

% mask mean beta
mb = double(mb);
mmb = mt * mean(mb(mb > 0));
mb(mb < mmb) = Inf;

% compute PSC
pscb = repmat(100 ./ mb, rma) .* double(b);

% remove illegal values
pscb(isinf(pscb) | isnan(pscb)) = 0;
